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abstract  IConllnua  on  tavataa  alda  II  naeaaaaff  and  IdattMly  by  bloek  inmbmtj 

A compact  shear  fracture  specimen  has  been  proposed  for  the  study  of  the  edge 
sliding  mode.  Mode  II,  of  crack  propagation.  The  purpose  of  this  investiga- 
tion was  to  numerically  determine  the  effect  that  crack  length,  specimen 
geometry,  and  applied  load,  has  on  the  Mode  II  stress  intensity  factors  for 
this  specimen.  Numerical  results  wore  compared  to  those  obtained  independent] 
by  another  investigator  using  boundary  collocation  and  photoelastic  techniquof 
In  addition,  the  initial  angle  of  propagation  was  numerically  determined  for- 
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selected  specimen  configurations  and  crapk  lengths.  The  Mode  II  stress 
intensity  factors  were  generated  by  using  displacement  and  strain 
energy  release  rate  methods.  The  initial  angle  of  propagation  was 
obtained  by  the  strain  energy  density  technique. 


The  optimum  specimen  configuration  established  dicing  this  investi- 
gation has  a tang  width  H = 1.0  inch.  This  mod^  configuration 
exhibited  a relatively  constant  nondimensionalized  Mode  II  stress 
intensity  factor  for  varying  crack  lengths.  For  a tang  width  of 
0.5  inch,  the  nondimensionalized  Mode  II  stress  intensity  factors 
were  found  to  be  load  dependent  for  short  crack  lengths.  For  a 
tang  width  H = 1.5  inch,  the  nondimensionalized  Mode  II  stress 
intensity  factors  were  Influenced  by  bending  of  the  specimen  over 
the  entire  range  of  crack  lengths  studied.  Boundary  effects  were 
also  noticed,  for  all  specimen  configurations  studied,  as  the  crack 
tip  approached  the  upper  and  lower  boundaries  of  the  specimen..  The 
numerically  determined  initial  angle  of  propagation  was  found  <to  be 
77°  for  the  specimen  configurations  and  crack  lengths  selected  be 
investigated.  This  study  also  showed  the  usefulness  of  the  finite 
element  method  in  determining  the  Mode  II  stress  intensity  factors 
for  different  loading  conditions  and  model  geometries. 
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CHAPTER  I 


INTRODUCTION  TO  PROBLEM 
1.1  General  Introduction 

The  use  of  high-strength  materials  in  the  design  of  engineering 
structures  has  lead  to  the  development  of  theories  which  can  be  used  to 
predict  the  reduced  strength  of  these  structures  caused  by  induced  or 
inherent  flaws  in  the  material.  Engineering  fracture  mechanics 
utilizes  the  concepts  of  stress  intensity  factor,  K , and  critical 
stress  intensity  factor,  , to  predict  this  reduced  strength.  The 
stress  intensity  factor,  K , is  a function  of  applied  load  and  geometry, 
while  the  critical  stress  intensity  factor,  , is  an  experimentally 
determined  constant  for  a given  material  and  mode  of  deformation. 

There  arc  three  possible  inodes  of  deformation  associated  with  a 
crack  as  shown  in  I’igure  1:  opening  mode.  Mode  1;  edge  sliding  mode, 

Mode  II;  and  tearing  mode.  Mode  III.  Until  recently.  Mode  I has  been 
considered  as  the  most  significant  mode  of  failure.  As  a conse<iuence, 
preceding  investigations  in  fracture  mechanics  have  dealt  primarily 
with  this  mode,  and  data  on  critical  stress  intensity  factors  is 
restricted  to  Mode  I loading.  Recent  investigations,  however,  indicate 
tliat  Mode  XI  may  be  a significant  mode  of  failure  in  certain  cases. 

Jones  .and  Chisholm  [1]  established  a compact  shear  fracture  specimen  to 
study  the  phenomenon  of  Mode  II  fracture  and  determine  accurate  critical 


stress  intensity  data. 
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The  object  of  this  investigation  was  to  check  the  performance  of 
the  compact  shear  fracture  specimen  by  numerically  generating  stress 
intensity  factors  using  finite  element  techniques.  In  addition,  initial 
crack  propagation  angles  were  numerically  predicted, 

1.2  Purpose  of  Investigation 

This  study  was  conducted  to  investigate,  by  finite  element  tech- 
niques, the  behavior  of  the  compact  shear  fracture  specimen  developed 
by  Jones  and  Chisholm,  In  particular.  Mode  II  stress  intensity  factors 
and  stress  boundary  conditions  were  numerically  generated  fur  comparison 
with  those  obtained  by  Jones  using  boundary  collocation  and  photoelastic 
methods.  In  addition,  the  initial  angle  of  propagation  which  could  not 
be  obtained  by  Jones  was  determined  for  selected  specimen  geometry  and 
crack  lengths, 

1.3  Scope  of  InvestiRation 

This  study  utilized  the  compact  shear  fracture  specimen 
established  by  Jones  and  Chisholm  to  numerically  determine  the  effect 
of  changing  crack  length,  applied  load,  and  specimen  geometry  on  the 
Mode  II  stress  intensity  factor,  . The  specimen, shown  in  Figure  2, 

has  a specluien  height  W , thickness  B , crack  length  a and  tang 
width  H . The  model  parameters  varied  in  this  investigation  were  crack 
length  a , and  tang  width  H . 

Point  and  uniform  loads  were  utilized  in  this  investigation  to 
simulate  conditions  that  could  be  achieved  in  a rigid  loading  frame, 
and  also  to  approximate  as  closely  as  possible  the  loading  conditiotts 
adopted  by  Jones.  The  stress  distribution  along  the  upper  boundary  and 


Figure  2.  Compact  Shear  Fracture  Specimen 


vertical  axis  of  symmetry  obtained  for  selected  load  conditions  and 
model  configurations  were  compared  to  those  assumed  in  the  boundary 
collocation  study. 
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The  initial  angle  of  propagation  of  the  crack,  for  a given  crack 
length,  applied  load,  and  model  geometry,  was  also  determined. 


1.4  Approach  Used  in  the  Investigation 


A numerical  solution,  the  finite  element  method,  was  employed  to 
generate  displacements  along  the  crack  flanks,  the  normal  and  tangential 
stress  distribution  along  the  upper  tang  boundary,  and  the  normal  stress 
distribution  along  the  vertical  axis  of  symmetry.  A plane  strain  con- 
dition was  utilized  for  the  numerical  analysis  throughout  this 
investigation. 

The  stress  intensity  factors  for  the  various  configurations 
studied  were  determined  by  the  displacement  [2]  and  strain  energy 
release  rate  methods  [2].  Numerical  results  were  compared  to  those 
obtained  in  a closed  form  solution  for  a model  with  similar  configura- 
tion and  loading  conditions,  and  to  the  bowidary  collocation  results. 

A dlscrlption  of  the  closed  form  solution  is  given  in  Appendix  H. 

The  numerically  generated  normal  and  tangential  stresses  along 
the  upper  boundary  and  vertical  axis  of  symmetry  represented  the  stress 
boundary  conditions  actually  occurring  in  the  model  for  a given  load 
and  model  configuration.  These  were  compared  to  the  stress  boundary 
conditions  assumed  by  Jones  and  Chisholm. 

A strain  energy  density  technique  (6]  was  used  to  numerically 
obtain  the  initial  angle  of  propagation  of  the  crack  for  selected  crack 
lengths  and  specimen  geometry. 


CHAPTER  II 


I ‘ 


FINITE  ELEl-tENT  TECHNIQUES  APPLIED  TO  LINEAR 
ELASTIC  FRACTURE  MECHANICS 

2.1  Introduction 

The  finite  element  method  was  utilized  to  generate  displacements 
along  the  free  flanks  of  the  crack,  and  the  normal  and  tangential  stress 
distributions  along  the  vertical  axis  of  syuanetry  and  upper  boundary  of 
the  compact  shear  fracture  specimen  under  consideration.  The  methods 
used  to  calculate  Mode  H stress  intensity  factors  for  given  specimen 
configuration,  crack  length  and  applied  loading  were  the  displacement 
methc'd  and  the  method  of  strain  energy  release  rate.  A brief  descrip- 
tion of  the  Unite  element  program  used  is  given  in  Appendix  A.  The 
initial  angle  of  propagation  for  selected  cases  was  numerically  determined 
by  the  strain  energy  density  method. 

2.2  Displacement  Method 

The  displacement  method  [2]  utilized  the  displacements  of  nodal 
points  arong  the  crack  flanks  and  equations  describing  the  displacement 
field  near  the  crack  tip.  The  equations  describing  the  displacement 
field  are  those  derived  by  Uestetgaard  [3],  and  shown  below  for  Mode  I 
and  Mode  II  plane  strain  deformations.  For  Mode  I loading, 
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while  for  Mode  II  loading. 


y ^L~ 


2v  + cos 


0 ] " 
2J  _ 


^11  r^T 

U 2tt 


-1  + 2v  + sin 


where  u is  the  shear  modulus,  V is  Poisson's  ratio,  r , 0 , u 

and  V are  defined  in  Figure  3.  Kj  and  are  the  Mode  1 and 

Mode  II  stress  Intensity  factors. 

By  using  nodal  displacements,  u*  and  v*  , along  the  crack 
flank  as  determined  by  the  finite  element  method,  a stress  Intensity 

'h  ^ 

factor  or  can  be  found  at  each  nodal  point  by  use  of 

Equations  (2.1)  through  (2.4).  Chan  [2]  found  that  the  most  accurate 
* * 

values  of  and  are  attained  by  using  the  equations  for 

and  u^^  with  6 « 180''  . Thus, 


4(1  - v^) 


(2Tt)^^^E  1 


4(1  - v‘) 


1/2  • 
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If  the  exact  displacements  were  used  as  r approaches  zero,  the  exact 
stress  intensity  will  be  determined.  Since  the  finite  element  method 
gives  accurate  solutions  along  the  flanks,  and  very  inaccurate  solutions 
at  the  crack  tip,  a tangent  extrapolation  must  be  used  to  determine 
stress  intensity  factors. 

This  is  accomplished  by  plotting  a graph  of  the  stress  intensity 

'ff 

factor,  K , versus  the  nondimensionalized  distance  from  the  crack  tip, 
r/a  , Figure  4,  using  Equations  (2.5)  and  (2.6).  A tangent  to  the 
straight  line  region  of  this  curve  is  then  extrapolated  back  to  the 
point  where  it  intersects  the  stress  intensity  axis.  This  intercept 
is  taken  as  the  value  of  the  stress  intensity  factor  K , A least 
squares  fit  to  the  data  in  the  straight  line  region  was  used  to 
determine  the  intercept  in  this  study.  As  can  be  seen  in  Figure  4,  the 
straight  line  region  does  not  extend  the  length  of  the  crack  flank. 

The  observed  nonlinearity  is  due  to  the  boundary  condition  imposed  at 
the  load  points,  and  the  fact  that  the  equations  for  displacement, 
Equations  (2.1)  through  (2.4),  are  strictly  valid  near  the  crack  tip. 

The  procedure  used  in  this  investigation  was  to  perform  a least  squares 
fit  in  the  region  of  r/a  ratios  ranging  from  0.1  to  0.2. 

2.3  Strain  Energy  Release  Rate 

The  concept  of  strain  energy  release  rate  [2]  states  that 
whenever  the  strain  energy  released  by  the  structure  is  greater  than 
the  energy  needed  to  create  new  crack  surface  area,  the  crack  will 
propagate  unstably.  Mathematically,  strain  energy  released  rate,  G , 


can  be  written  as; 


90 


I31NI  SSatJXS 
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G = ^ 

^ dA  » 


(2.7) 


where  U is  the  strain  energy  stored  in  the  structure  and  A is  the 
new  crack  surface  area  formed  as  the  crack  extends.  Strain  energy 
release  rate  is  related  to  Mode  I and  Mode  II  stress  intensity  factors 
by  the  relation  [4] 


2 2 


(2.8) 


1 - V 


This  equation  is  valid  for  plane  strain  conditions. 

The  finite  element  method  can  be  used  to  solve  for  the  strain 
energy  release  rate,  expressed  by  Equation  (2.7),  by  the  following  means. 
The  total  strain  energy  of  the  structure,  U^,  , is  determined  by 
numerically  summing  the  strain  energy  of  each  element  of  the  structure. 
The  crack  is  then  extended  by  a small  Incremental  length,  Aa  , and  the 
total  strain  energy,  U,^  * is  again  dotemined.  For  plane  strain 
analysis,  with  the  thickness  of  the  model  set  equal  to  unity.  Equation 
(2.7)  can  now  be  rewritten  as: 


u;  - 


(2.9) 


Utilizing  Equation  (2.9)  and  the  fact  that  for  the  fracture  specimen 
being  studied,  Mode  I stress  intensity  factors  were  negative  and 
therefore  could  be  neglected,  can  be  determined  from  the  relation: 


2 2 


<1  * 
1 - V 


Ui  - u^, 


(2.10) 


2.4 


Strain  Energy  Density 

Because  cracks  which  are  not  oriented  perpendicular  to  the  applied 
load  tend  to  propagate  in  a direction  other  than  along  the  axis  of  the 
crack,  Sih  [5]  proposed  the  concept  of  strain  energy  density  in  order  to 
analytically  determine  the  direction  of  crack  growth.  Sih  determined 
that  the  magnitude  of  the  energy  field  in  the  vicinity  of  the  crack  tip 
can  be  written  as: 


S = a^j^K^  + + a22Kjj 


(2,11) 


where 


and 


®11  “ Idil  * 


&12  = ® " (k  “ 1)] 


^22  ° ((k  + 1)(1  - cos  0)  + (1  + cos  0)(3  cos  0-1)]  , 

(2.12) 


In  these  equations,  k ■ 3 - 4v  for  plane  strain  and  (3  - \>)(1  + v) 
for  plane  stress  problems,  0 is  defined  In  Figure  3,  u is  the  shear 
modulus  and  and  are  the  Mode  I and  Mode  II  stress  Intensity 

factors.  The  concept  of  strain  energy  density  states  that  a crack  will 
propagate  in  the  direction  for  which  the  strain  energy  density  S 
possesses  a stationary  minimum  value,  or 
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H = 0 . (2.13) 

The  initial  direction  of  propagation  can  now  be  determined  by  first 
using  finite  element  techniques  to  solve  for  the  stress  intensity 
factors,  Kj,  and  . Equations  (2.11)  and  (2.12)  are  then  used  to 

numerically  identify  the  value  of  6 which  will  result  in  a minimum 
value  of  S . In  this  study,  the  solution  of  this  problem  produced  two 
values  of  0 which  yielded  a minimum  value  of  S with  the  correct 
value  being  determined  by  a physical  argument. 


CHAPTER  III 


MODEL  USED  IN  THE  INVESTIGATION  AND 
VERIFICATION  OF  NUMERICAL  PROCEDURE 

3.1  Introduction 

As  stated  earlier,  the  model  considered  in  this  study  Is  a 
fracture  specimen  proposed  by  Jones  and  Chisholm  and  previously  shown 
in  Figure  2.  This  model  was  gridded  for  generation  of  displacements 
and  stresses  by  the  finite  element  method.  The  boundary  and  load 
conditions  applied  to  the  model  are  discussed  in  the  following  para- 
graphs and  verification  of  the  gridding  program  and  finite  element 
technique  is  established. 

3.2  Grid  Pattern 

A typical  grid  pattern  utilized  ia  this  investigation  is  shown 

in  Figure  3.  The  pattern  consists  of  a fine  region  at  the  crack  tip, 

a transition  region  and  a course  region.  In  the  fine  region,  a typical 

2 

ratio  of  clement  area  to  the  square  of  crack  length,  A/a  , is 

_3 

2.6  X 10  for  a/U  ■>0.5  . This  ratio  for  the  course  region  is 
-2 

I.OA  X 10  . Ratios  were  chosen  in  each  region  to  assure  optimal 

convergence  to  the  true  solution  [ 2 ] . The  material  above  the  loading 
pins  was  not  included  in  the  finite  element  model  because  it  did  not 
contribute  any  significant  stiffness.  The  fine  grid  region  was  moved 
with  the  crack  tip  as  crack  length  was  varied.  This  required  that  the 


Figuro  5.  Typical  Grid  Fattern,  a/W  » 0.5 
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mesh  be  regenerated  each  time  a new  crack  length  was  considered.  In 
order  to  decrease  the  amount  of  time  spent  in  regridding  the  model,  a 
grid  generation  code  obtained  from  another  finite  element  program  (SAAS- 
Stress  Analysis  of  Axisymmetric  Solids) [6]  was  modified  to  grid  the 
specimen.  The  modification  and  required  irput  into  this  program  is 
discussed  in  detail  in  Appendix  C.  Briefly,  this  program  interpolates 
the  position  of  internal  nodal  points  from  those  specified  on  the  ex- 
ternal boundaries.  The  program  then  eliminates  nodaJ  points  as 
described  in  Appendix  C in  orde : to  achieve  the  desired  reduction  of  the 
grid  pattern.  Finally,  the  elerient  Information  is  obtained  by  using 
generated  nodal  point  data. 

3,3  Model  Geometry 

A list  of  the  parameters  considered  in  this  investigation  is 
given  in  Table  I.  The  ratio  of  crack  length  to  specimen  height,  a/U  , 
was  varied  from  0.1  to  0.8.  This  range  was  selected  to  determine  the 
effect  that  the  model  boundaries  have  on  the  Mode  II  stress  intensity 
factors  for  short  and  long  cracks.  The  dependence  of  the  Mode  II  stress 
intensity  factor  on  distance  from  the  load  point  to  crack  flank  was 
studied  by  changing  the  tang  width,  H , of  the  specimen.  Three  widths 
were  considered;  11  <■  l.S  inches,  H 1.0  inch,  and  H « 0.5  inch. 

For  H ■■  l.S  inches,  the  same  grid  pattern  utilised  for  II  “ 1.0  inch 
was  adopted  with  the  outer  boundary  elements  expanded  to  give  the  desired 
tang  width.  The  grid  pattern  for  H ■ 0.5  inch  had  to  be  regenerated. 
This  was  easily  accomplished  by  a toodification  of  the  griding  program 
which  allowed  it  to  accept  data  for  H ■>  1.0  inch,  and  grid  the  model 
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for  H = 0.5  inch.  The  modification  restricted  the  range  over  which 
the  program  will  search  for  nodal  points,  thus  giving  the  desired  tang 
width. 

3.4  Load  and  Boundary  Conditions 

The  loading  conditions  used  in  this  investigation  consisted  of 
point  and  uniform  loads  applied  to  the  specimen  as  shown  in  Figures  6a 
and  6b,  respectively.  In  order  to  simulate  the  loading  conditions  used 
in  the  boundary  collocation  analysis  [1],  it  was  necessary  to  pin  the 
load  points  in  the  y-direction  as  shown  in  Figure  6a.  The  left  boundary 
of  the  model  was  pinned  in  the  y-direction  because  it  is  a line  of 
symmetry.  The  uniform  loading  condition,  with  the  nodal  points  on  the 
upper  boundary  pinned  in  the  y-direction  as  shown  in  Figure  6b,  was 
selected  to  simulate  a fixed  grip  loading  frame.  The  crack  tip  was 
pinned  in  the  z-direction  to  eliminate  any  rigid-body  movement  of  the 
crack  flanks  in  the  z-direction. 

3 . 5 Verification  of  Griddlng  Program  and  Numerical  Techniques 

Verification  of  the  gridding  program  and  numerical  techniques 

for  obtaining  stress  intensity  factors  was  accomplished  by  generating 
a representative  grid  pattern  and  loading  the  compact  shear  specimen  as 
a compact  tension  specimen.  The  transformation  from  compact  shear  to 
compact  tension  specimen  was  easily  made  by  simply  applying  the  loads 
at  the  side  and  vertical  axis  of  symmetry  along  a line  perpendicular  to 
the  crack  axis  as  shown  in  Figure  7,  The  Mode  I stress  intensity  factor 
was  then  numerically  obtained  usiv g both  the  displacement  and  strain 
energy  release  rate  techniques  and  compared  to  an  available  closed  form 


Boundary 
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solution.  A closed  form  solution  for  the  compact  tension  specimen, 
shovm  in  Figure  7,  is  given  by  Sih  [10]  as: 


4 “ 


P(2W  + a) 


/rr  h(W  - a) 


3/2 


5.  ^ £ 

[W’  W’  bj 


(3.1) 


where  F | is  given  for  various  parameters  by  Sih  [10]  and  will 

be  used  here  only  for  the  geometry  considered.  The  boundary  conditions 
necessary  for  numerical  generation  of  the  Mode  I stress  intensity  factor 
are  shown  in  Figure  7.  '.’ith  these  boundary  conditions  and  the  model 
configuration  for  a compact  shear  fracture  specimen  with  a/W  =■  0.5  and 
H = 1.0  inch  , the  physical  parameters  of  Equation  (3.1)  become 
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A r-  ratio  of  0.3  was  chosen  because  this  was  the  closest  value  listed 
b 

by  Sih  which  fit  the  loading  conditions  applied  to  the  compact  tension 
specimen.  It  also  should  be  noted  that  the  value  of  stress  intensity 
factor  given  by  Equation  (3.1)  is  not  the  same  as  that  determined  by 
the  displacement  tectuiique.  The  relation  between  these  two  values  is 
given  by 


Kj  “ /it  kj. 


I 


(3.2) 
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where  is  determined  by  Equation  (3.1)  and  is  determined  by  the 
displacement  method. 

A graph  of  stress  intensity  factor  versus  nondimensionalized 
distance  from  the  crack  tip  for  the  compact  tension  specimen  is  given 
in  Figure  8.  A least  squares  fit  was  performed  on  the  points  in  the 
region  r/a  from  0.1  to  0.2  and  , the  Intercept,  was  found  to  be 
Kj  = 500  psi  /in  . The  value  determined  from  Equations  (3.1)  and  (3.2) 


with  F 


£ ^ JC 

W’  W’  b 


= 2.010  is  = 493  psi  /Tn  . These  values  are 
in  good  agreement  with  each  other  and  therefore  verify  the  displacement 
technique  and  grid  generation  program. 


CHAPTER  IV 


ANALYSIS  OF  RESULTS 


4.1  Introduction 

Having  verified  the  numerical  methods  and  grid  generation  program 
used  in  this  investigation,  the  behavior  of  the  compact  shear  fracture 
specimen  as  a function  of  crack  length,  loading  condition  and  specimen 
geometry  is  discussed.  Mode  II  stress  intensity  factors  generated  by 
finite  element  techniques  are  compared  to  those  determined  by  boundary 
collocation  analysis,  and  to  a closed  form  solution  for  geometry  and 
loading  conditions  that  closely  approximate  those  of  the  compact  shear 
fracture  specimen.  The  numerically  generated  stress  distribution  on  the 
upper  model  boundary  and  vertical  axis  of  symmetry  are  compared  with  those 
assumed  by  Jones  and  Chisholm  [1]  in  a boundary  cfllocatlon  study.  The 
initial  angle  of  propagation  is  also  given  for  sel acted  crack  lengths 
and  specimen  tang  widths. 

4.2  Stress  Intensity  Factors 

Numerically  generated  Mode  II  stress  intensity  factors  for 
selected  nondimensionalized  crack  lengths  and  tang  widths  are  listed  in 
Table  II,  This  data  is  graphically  portrayed  in  Figure  9 where  the 
Mode  II  stress  intensity  factors  have  been  nondimonsionalized  by 
dividing  by  the  normal  stress  in  the  tang,  0 » P/BH  , and  the 

square  root  of  crack  length,  /a  , The  curves  shown  in  Figure  9 
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TABLE  II 

NUMERICALLY  GENERATED  MODE  II 
STRESS  INTENSITY  FACTORS 

Mode  II  Stress  Intensity  Factors  (psi  / In . ) 
a/W  H = 0.5  Inch  H = 1.00  Inch  H = 1.50  Inches 

Displacement  Energy  Displacement  Energy  Displacement  Energy 

0.1  - - 63.29  58.35 

0,2  148.09  137.30  78,59  69,00  97.21  46.52 

0.3  152.49  145.10  92.80  92.50 

0.4  - - 115.83  107.00  101.69  92.90 

0.5  137.68  136.70  122.06  118.00  114.71  107.00 

0.6  - - 145.85  139.00  142.41  129.50 

0.7  158.44  166.44  173.09  164.00  169.48  152.40 

0.8  - - 209.27  199.79 


a/W  Ration  for  Selected 


27 


indicate  that  the  specimen  having  a tang  width  H = 1.0  inch  is  the 
most  stable,  as  evidenced  by  a relatively  constant  nondimens ionalized 
Mode  II  stress  intensity  factor.  The  curves  for  H = 0.5  inch  and 
H = 1.5  inches  exhibit  rather  wide  variations  in  the  nondimensionalized 
Mode  II  stress  intensity  factors  with  varying  crack  length.  These  cases 
will,  therefore  be  examined  in  more  detail. 

Figure  10  shows  the  variation  of  the  nondimensionalized  Mode  II 
stress  intensity  factor  with  crack  length  for  a tang  width  H « 0.5  inch. 
Excellent  correlation  is  shown  between  results  obtained  by  finite 
element  techniques  and  those  determined  from  a closed  form  solution  [7]. 
However,  the  boundary  collocation  data  does  not  agree  with  the  numerical 
and  closed  form  results.  The  curves  for  the  finite  element  and  closed 
form  solutions  show  a region  of  constant  nondimensionalized  Mode  II 
stress  intensity  factors  from  a/W  = 0.4  to  a/W  =0.7  . For  a/W 
ratios  greater  than  0.7,  the  lower  specimen  boundary  influencies  the 
Mode  II  stress  intensity  factors.  At  a/W  ratios  lower  than  0.4,  the 
curves  show  a load  dependence  as  the  load  is  now  being  applied  close  to 
the  crack  tip.  The  boundary  collocation  curve  exhibits  constant  non- 
dimensionalized Mode  II  stress  Intensity  factors  over  the  entire  range 
of  a/W  ratios  because  of  incorrect  stress  boundary  conditions  assumed 
in  the  boundary  collocation  analysis.  This  error  in  assumed  stress 
boundary  conditions  will  be  discussed  in  detail  in  a later  section. 

Figure  11  shows  the  variation  of  nondimensionalized  Mode  II 
stress  intensity  factors  with  crack  length  for  a tang  width  H = 1.0 
inch.  All  curves  are  in  good  agreement  with  each  other  and  have  a 
region  of  constant  nondimensionalized  Mode  II  stress  Intensity  factors 


Nondiraeosionallzed  Mode  II  Stress  Intensity  Factor 
Versus  a/W  Ratio  for  a Tang  Width  K = 0.5  Inches 
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from  a/W  = 0.2  to  a/W  = 0.6  . At  a/W  ratios  above  0.6  and  below 
0.2,  the  specimen  boundaries  are  being  sensed  by  the  crack  tip,  causing 
the  curves  to  turn  upward. 

Figure  12  gives  the  variation  of  nondimensionalized  Mode  II  stress 
intensity  factors  with  crack  length  for  a tang  width  H = 1.5  inches. 

The  finite  element  and  closed  form  results  do  not  exhibit  a region  over 
which  the  nondimensionalized  Mode  II  stress  Intensity  factors  are 
constant  because  the  specimen  is  now  behaving  like  a beam  in  three  point 
banding.  The  influence  of  bending  is  evident  at  all  a/W  ratios 
particularly  chose  less  than  a/W  = 0.3  , for  in  this  region  the  curves 
were  expected  to  turn  upward  as  the  crack  tip  sensed  the  upper  boundary 
of  the  specimen.  The  boundary  collocation  data  exhibits  a different 
behavior  from  that  of  the  finite  element  and  closed  form  results  because 
of  incorrect  stress  boundary  conditions  assumed  in  the  boundary  colloca- 
tion study. 

4.3  Stress  Boundary  Conditions 

Stress  boundary  conditions  were  determined  along  the  upper 
boundaries  of  each  tang  and  along  Che  vertical  axis  of  symmetry  in  order 
to  check  those  assumed  in  the  bo^indary  collocation  analysis  and  to 
provide  useful  data  for  future  investigations  utilizing  boundary  value 
techniques.  Of  particular  interest  in  this  investigation  was  the  stress 
distribution  along  Che  upper  tang  boundary  and  vertical  axis  of  symnetry 
for  application  of  a point  load  with  the  load  points  fixed  in  a direction 
perpendicular  to  Che  direction  of  loading,  as  shown  in  Figure  6a.  This 
method  of  loading  was  of  interest  because  it  simulates  that  used  by 


Nondlniensional i zed  Kode  II  Stress  Intensity  Factors 
Versus  a/W  Ratio  for  a Tang  Width  H = 1.5  Inches 
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Jones  and  Chisholm  so  that  resultant  numerically  generated  stress  dis- 
tributions could  be  used  to  check  the  stress  boundary  conditions  assumed 
by  Jones  and  Chisholm.  A uniform  load  with  each  nodal  point  on  the 
upper  tang  boundary  fixed  in  a direction  perpendicular  to  the  direction 
of  loading,  as  shown  in  Figure  6b,  was  also  applied  in  order  to  determine 
the  stress  boundary  conditions  needed  in  order  to  perform  boundary  value 
analyses  for  this  common  type  of  loading.  The  two  loading  conditions 
just  discussed  are  obtainable  through  a rigid  loading  frame.  The  stress 
distribution  for  two  other  loading  conditions  are  also  discussed  in  this 
section.  These  loading  conditions  are  a point  load  with  all  nodal  points 
on  the  upper  tang  boundary  fixed  in  a direction  perpendicular  to  the 
direction  of  loading,  and  a uniform  load  with  the  center  point  of  each 
tang  fixed  in  a direction  perpendicular  to  the  direction  of  loading. 
Although  these  conditions  are  easily  simulated  in  a numerical  study, 
they  are  very  difficult  to  obtain  experimentally  and  were  considered  in 
this  investigation  simply  to  determine  the  sensitivity  of  stress 
intensity  factors  to  changes  in  boundary  conditions. 

The  stress  boundary  conditions  assumed  in  the  boundary  collocation 
analysis  are  shown  in  Figure  13.  The  stress  distribution  along  the  upper 
tang  boundary  consists  of  a uniform  tensile  stress  on  the  outside  tang, 
a uniform  compressive  stress  a^ong  the  center  tang  and  a cosine  distribu- 
tion of  tangential  stress  along  the  top  of  both  tangs.  A bilinear 
bending  stress  was  assumed  along  the  vertical  axis  of  symmetry. 

The  numerically  generated  stress  distributions  along  the  upper 
tang  boundary  and  vertical  axis  of  symmetry  for  a/W  •=  0.5  , H “ 1.0 
inch  and  plane  strain  conditions  are  shown  in  Figures  14  through  19  for 
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Nus^e^ically  Generated  Stress  Distribution  Along  Upper 
Tang  Boundary  for  a/W  = 0.5,  H = 1.0  Inch,  Point  Load 
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selected  load  and  displacement  boundary  conditions.  Figures  14  through 
16  show  the  stress  distributions  for  a point  and  uniform  load  with  the 
center  point  of  each  tang  fixed  in  a direction  perpendicular  to  the 
direction  of  loading.  Figures  17  through  19  show  the  stress  distribu- 
tions for  a point  and  uniform  load  with  each  nodal  point  along  the  upper 
tang  boundary  fixed  in  a direction  perpendicular  to  the  direction  of 
applied  load. 

As  can  be  seen  by  referring  to  Figures  14  and  17,  the  normal  and 
tangential  stress  distribution  along  the  upper  tang  boundaries  caused 
by  the  application  of  a point  load  does  not  change  as  the  displacement 
boundary  conditions  are  varied.  Figures  15  and  18  show  that  the 
tangential  stress  distribution  along  the  upper  tang  boundaries  caused 
by  the  application  of  a uniform  load  is  in  fact  the  only  distribution 
affected  by  variation  in  displacement  boundary  conditions.  Figures  16 
and  19  show  practically  no  change  in  the  normal  stress  distribution  along 
the  vertical  axis  of  symmetry  for  variation  in  displacement  boundary 
conditions. 

Comparison  of  Figures  13,  14  and  17  Indicates  that  the  numerically 
generated  stress  distributions  for  point  loading  are  not  in  good 
agreement  with  those  assumed  in  the  boundary  collocation  analysis.  The 
major  differences  occur  in  the  normal  stress  distribution  along  the 
upper  tang  boundary.  This  difference  in  stress  distribution  is  respon- 
sible for  the  discrepancies  observed  in  the  nondimens ionallzed  Mode  II 
stress  Intensity  factors  determined  by  numerical  techniques  and  boundary 
collocation  analysis. 
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The  stress  distributions  caused  by  the  application  of  a uniform 
load  with  all  the  nodal  points  along  the  upper  tang  boundary  fixed,  as 
shown  in  Figure  18,  were  somewhat  similar  to  those  assumed  in  the 
boundary  collocation  analysis,  shown  in  Figure  13.  The  major  difference 
in  the  two  distributions  is  that  the  numerically  generated  tangential 
stress  distribution  on  the  center  tang  is  not  a cosine  distribution  as 
assumed  in  the  boundary  collocation  analysis.  Figures  16  and  19  show 
that  the  bending  stress  along  the  vertical  axis  of  symmetry  for  both 
point  and  uniform  loading  is  not  a simple  bilinear  distribution  as  shown 
in  Figure  13. 

Numerically  generated  stress  distributions  were  also  obtained  for 
other  a/W  ratios  and  tang  widths  resulting  in  distributions  having  the 
same  form  but  different  magnitudes  then  those  shown  in  this  section. 
Plane  stress  conditions  were  also  tried  for  selected  load  conditions  and 
model  geometries  resulting  in  no  difference  in  stress  distributions  when 
compared  to  those  determined  using  plane  strain  cor'^itions.  The  Mode  II 
stress  intensity  factors  were  unaffected  by  changes  in  loading  and 
displacement  boundary  conditions  discussed  in  this  section. 

4.4  Initial  Angle  of  Crack  Propagation 

The  angle  of  propagation  6 is  defined  in  Figure  3.  By  use  of 
the  strain  energy  density  technique  proposed  by  Slh  [5]  and  discussed 
in  Chapter  II,  the  angle  of  propagation  was  determined  for  a/W  >0.5 
through  0.8  with  H =>  1.0  inch  . These  values  are  listed  in  Table  III. 
These  particular  values  of  crack  length  and  specimen  tang  width  were 
chosen  because  the  experimental  work  performed  by  Jones  considered  an 


TABLE  III 


NUMERICALLY  PREDICTED  INITIAL  ANGLE  OF 
PROPAGATION  FOR  SELECTED  a/W  RATIOS  WITH 
A TANG  WIDTH  H = 1.00  INCH 


a/W 

THETA  0 

0.5 

77“ 

0.6 

ir 

0.7 

77“ 

0.8 
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a/W  .ratio  of  0.8  with  specimen  tang  width,  H = 1.0  inch  . With  these 
parameters,  and  notching  the  model  along  the  plane  of  the  crack,  Jones 
and  Chisholm  were  able  to  obtain  a straight  fracture  of  the  specimen, 
0=0°.  As  can  be  seen  from  Table  III,  the  strain  energy  density 
technique  predicts  that  the  initial  angle  of  propagation  is  75°  for 
a/W  = 0.8  and  77°  for  a/W  = 0.5  to  0.7  . The  path  of  the  crack  was 
not  investigated  after  the  initial  angle  of  propagation  because  of  the 
amount  of  computer  time  necessary  to  accomplish  this. 


CHAPTER  V 


CONCLUSIONS  AND  RECOMMENDATIONS 


5.1  Conclusions 

The  effect  of  crack  length  and  specimen  geometry  on  the  Mode  II 
stress  intensity  factors  for  the  compact  shear  fracture  specimen  shown 
in  Figure  2 was  investigated.  Stress  intensity  factors  were  generated 
by  the  displacement  and  strain  energy  release  rate  methods  and  compared 
to  those  determined  by  boundary  collocation.  Both  numerical  and 
boundary  collocation  results  obtained  for  a tang  width  H = 1,0  inch 
show  that  nondimensionalized  Mode  II  stress  intensity  factors  are 
constant  for  this  specimen  over  a wide  range  of  crack  lengths.  This 
configuration  of  the  compact  shear  fracture  specimen  will  therefore  be 
most  suitable  for  determining  critical  Mode  II  stress  intensity  factors, 

^Ic  • 

Numerical  results  for  U 0.5  inch  and  H « 1.5  inches  show 
that  nondimensionalized  Mode  II  stress  intensity  factors  are  not  as 
stable  with  Increasing  crack  length  as  the  boundary  collocation  results 
suggest.  The  stress  intensity  factors  obtained  for  a tang  width 
H 1.5  inches  are  affected  by  bending  of  the  specimen,  a fact  that 
was  not  brought  out  by  the  L'^undary  collocation  results.  Results  for 
tang  width  U » 0.3  inch  indicated  that  the  nondimensionalized  Mode  II 
stress  Intensity  factors  were  influenced  by  applied  load  for  a/W  ratios 
less  than  0.4  . Because  of  these  nonlinearities,  the  conflguratlo..s 


45 


having  tang  width  H = 0.5  inch  and  H = 1.5  inches  are  not  suitable 
ior  determining  the  critical  Mode  II  stress  intensity  factor,  • 

For  the  crack  lengths  and  specimen  geometries  investigated,  the  non- 
dimensionalized  Mode  II  stres  intensity  factors  were  found  to  be 
independent  of  the  applied  load,  and  uninfluenced  by  the  lower  free 
boundary  through  the  intermediate  ranges  of  crack  lengths. 

Numerical  results  indicated  that  the  tangential  stress  along  the 
upper  boundary  of  the  center  tang  did  not  follow  a cosine  distribution 
as  assumed  in  the  boundary  collocation  analysis.  In  addition  the 
numerically  determined  bending  stress  distribution  along  the  vertical 
axis  of  symmetry  was  not  found  to  be  the  simple  bilinear  distribution 
which  was  used  in  the  boundary  collocation  investigation.  Therefore, 
incorrect  stress  boundary  conditions  assumed  in  the  boundary  collocation 
analysts  were  responsible  for  differences  observed  in  the  behavior  of 
Che  ncridimensionaltzed  Mode  II  stress  intensity  factors  when  numerical 
and  collocation  results  are  compared. 

Iho  initial  angle  of  propagation  was  found  numerically  to  be  11" 
for  rhe  geometry  and  crack  lengths  considered.  This  differs  from  the 
experimental  wolk  performed  by  Jones  and  Chisholm  (1),  which  showed  an 
angle  of  O'  of  crack  propagation.  The  snglc  of  0'  might  have  been  the 
result  ct  notching  the  fracture  specimen  along  the  plane  of  the  crack. 
The  difference  between  the  numerically  predicted  angle  of  propagation 
and  experimental  "esult-S  might  not  have  been  so  groat  if  it  wore 
possible  to  numerically  determine  the  complete  trajectory  of  the  crack. 
This  was  not  done  because  of  the  amount  of  computer  time  involved. 
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This  study  also  showed  the  usefullness  of  the  finite  element 
method  in  determining  Mode  II  stress  intensity  factors,  for  different 
loading  conditions  and  specimen  parameters.  The  accuracy  of  the  finite 
element  techniques  discussed  in  this  investigation  was  brought  out  by 
the  good  agreement  with  a closed  form  solution. 

5.2  Recommendations 

It  is  recommended  that: 

1)  An  experimental  verification  of  the  stress  boundary  conditions 
and  Mode  II  stress  intensity  factors  predicted  by  the  finite  element 
techniques  used  in  this  investigation  bi  performed. 

2)  A verification  of  the  crack  trajectory  by  experimental  or 
numerical  techniques,  with  a notch  Introduced  into  the  model  for 
comparison  with  the  experimental  work  performed  by  Jones  and  Chisholm 


be  undertaken. 
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APPENDIX  A 

FINITE  ELEMENT  PROGRAM  USED  IN  THE  INVESTIGATION 
A.l  Finite  Element  ProRram 

The  finite  element  program  utilized  In  this  Investigation  was 
SSAP-2  (Static  Analysis  Program  For  Solid  Structures).  A complete 
listing  of  the  program  as  well  as  a description  of  the  elements  and 
Input  to  the  program  Is  given  In  the  program  manual  [8] . A brief 
description  of  the  plate  element  used  In  the  study  will  be  given  here. 

A.  2 Isoparametric  Elements 

A isoparametric  element  (9]  is  defined  as  one  where  both  the 
displacement  and  geometry  of  the  element  are  described  by  the  same 
parameter.  This  means  that  the  relation  between  the  local  and  global 
coordinate  systems  as  well  as  the  displacement  approximation  for  the 
element  are  given  in  terms  of  the  same  interpolation  function.  Air 
interpolation  function  is  one  that  has  a unit  value  at  one  nodal  {xiiut 
and  is  zero  at  all  others  in  the  element.  Two  advantages  of  using 
interpolation  functions  are: 

1)  If  continuity  of  geometry  and  displacement  both  within  and 
between  adjacent  elements  are  satisfied,  compatibility  is  satisfied  in 
global  coordinates. 

2)  If  the  interpolation  function  is  able  to  give  rigid  body 
displacements  in  the  local  coordinate  system,  rigid  body  displacement 
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and  constant  strain  are  satisfied  in  the  global  system.  This  is 
necessary  for  convergence  of  the  method. 

The  interpolation  functions  used  in  GSAP  will  now  be  discussed. 
For  the  quadrilateral  element  shown  in  Figures  20  and  21  the  local  and 
global  coordinates  are  related  by 


4 

X =-■  1 h.x 

i=l  ^ ^ 

and 

4 

• 

i=l 


(A.l) 


where  the  interpolation  functions, 

hj^  “ 

h2  = -J-  (1  + s)(l  - t) 

h^  =■  (1  + s)(l  + t) 

and  . 

(1  - s)(l  + t) 


h^  , are  given  as: 


$ 


f 

(A. 2) 


where  s and  t are  defined  In  Figure  21. 

The  displacements  of  the  clement  are  written  in  terms  of  the 
same  Interpolation  functions  as  follows: 


Ujj(s,t) 
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h.u  . 
1 xi 
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and 


Uy(s,t)  = 


4 

Z 

i=l 


h.u 
i yi 


The  element  strains  can  now  be  written  as 


and 
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These  equations  can  be  written  In  matrix  form  as: 
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where  and  ^ are  the  displacements  In  the  x and 
respectively  and 


and 


where 


and 


U - [h,  h,  h,  h,  ] 

,x  l,x  2,x  3,x  4,x' 
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directions. 
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Since  the  are  given  in  terms  of  the  natural  coordinates,  s and 
t , the  chain  rule  must  be  applied  in  order  to  compute  the  derivatives 
in  terms  of  the  x and  y coordinate  system.  The  result  of  the 
application  of  the  chain  rule  is  shown  in  matrix  form  below  as: 


at  y ^ ~y 

,x  ,x  ^ ,t  ^,8 


s t 

>y  .y. 


-X  X 

»t  , s 


(A.  7) 


where  J is  the  Jacobian  which  is  defined  by 


■ ^s ''.t  - *.t  '.S  • 


X = Z h.  X.  , 

,s  i,s  i 


' t " ^ ^i  t ’^i  » 

i«l  ^ 


(A.8) 


y » r.  h.  y . 

i.s  'i 


Dow  Equations  t^A.7)  and  (A.8)  can  be  used  with  Equations  (A. 6)  to 
determine  the  straia'displacement  matrix  of  Equation  (A. 5). 


Element  Stiffness 


The  element  stiffness  for  unit  thickness  is  given  as 


K “ a £ 

4 


a dA 


(A.9) 


jc  is  the  stress-strain  matrix  and  the  integration  is  performed  over  the 
area  of  the  element.  This  equation  can  be  written  in  terms  of  the 
natural  coordinate  system  as 
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rl  rl 
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c a J ds  dt 


(A.  10) 


Standard  numerical  integration  is  used  to  determine  K for  a given 
element . 


A. 4 Total  Stiffness  Matrix 

Once  the  element  stiffnesses  have  been  determined,  the  total 
stiffness  matrix  is  obtained  by  summing  element  stiffnesses  in  the 
conventional  manner. 


^OTAL  “ ^ Element 


(A. 11) 


A complete  description  of  the  finite  element  method  is  given  by  Deal  [9]. 
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APPENDIX  B 

CLOSED  FORM  SOLUTION 
B.l  Geometry  and  Equations 

The  verification  of  the  finite  element  results  was  accomplished 
by  comparison  with  a closed  form  solution  for  a model  with  geometry  and 
loading  conditions  similar  to  those  of  the  compact  shear  fracture 
specimen  being  studied.  Since  a geometry  and  loading  condition  that 
exactly  matched  this  model  was  not  available,  a configuration  was 
considered  that  would  give  a lower  bound  to  the  finite  element  results. 
The  loading,  geometry  and  closed  form  solution  were  chosen  from  a 
handbook  by  Tada  [7j  and  the  model  is  shown  in  Figure  22.  This 
configuration  was  chosen  because  the  loading  applied  is  similar  to  the 
loading  on  the  compact  shear  fracture  specimen.  The  configuration  will 
give  a lower  bound  to  the  problem  because  it  is  an  infinite  strip  and 
therefore  stiffer  than  the  actual  plate.  The  equations  for  the  Mode  XI 
stress  intensity  factor,  , are  given  below: 
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APPENDIX  C 


GRID  GENERATION 


C.l  Introdviction 

Tlie  mesh  generation  program  utilized  in  this  investigation  is  a 
modified  version  of  one  used  in  SAAS  [6].  The  modification  is  in  the 
element  generation  subroutine  and  nodal  point  interpolation.  Element 
generation  was  changed  in  order  to  generate  triangular  elements  in 
regions  of  transitioti  from  course  to  fine  gridding  and  nodal  point 
interpolation  was  modified  to  make  grid  reduction  possible. 

C.2  Input  Required 

Input  into  the  gridding  program  will  be  described  by  considering 
the  gi'iddlng  of  a plate  with  a notch  as  shown  in  Figure  23.  The 
requirements  for  gridding  a crack  problem,  without  the  use  of  special 
crack  tip  elements,  are  that  a fine  grid  be  placed  in  the  crack  tip 
region  to  accommodate  the  high  stress  gradient  there.  Also,  in  order 
to  save  computer  time  and  storage,  a relatively  course  grid  should  be 
applied  at  the  boundaries.  These  two  requirements  can  be  accomiJlished 
by  a transition  region  of  triangular  elements.  The  procedure  for 
obtaining  input  data  for  this  program  is  as  follows.  Determine  what 
type  of  grid  pattern  is  needed  at  the  crack  tip,  using  only  quadrat- 
lateral  elements.  This  will  give  a very  fine  element  distribution  as 
shown  in  Figure  24.  Each  horizontal  line  with  the  same  z-coordinate  is 


M 1=2  1=3  1=4 1=9  1=10 


I=ll  R 


Figure  24. 


Grid  Pattern  for  Notched  Plate  Wthout  Triangular 
Transition  Elements 
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assigned  a J index  to  identify  it,  starting  with  the  z = 0 line.  The 
same  is  done  with  every  vertical  line  with  the  same  R coordinate,  except 
these  are  assigned  I indices;  therefore,  every  intersection,  nodal  point, 
has  an  I,J  index  to  identify  it.  The  transition  can  now  be  accomplished 
by  eliminating  sections  of  J rows  and  I columns  in  order  to  achieve  the 
desired  grid  pattern.  The  only  restriction  is  that  a section  of 
eliminated  J row  cannot  intersect  a portion  of  I column  that  will  be 
eliminated  because  the  program  will  not  be  able  to  generate  a valid 
element.  Figure  25  shows  a completed  grid  pattern. 

The  sequencing  of  the  input  data  will  now  be  described.  A listing 
of  input  data  for  this  problem  is  given  in  Figure  26.  The  first  three 
cards  are  control  cards  and  are  described  in  the  manual  for  SAAS  and 
the  control  data  for  the  finite  element  program  has  been  eliminated. 

The  cards  and  formating  are  described  below. 

Title  Card 

Describe  grid  pattern 

Format  (8A10) 

Columns  1-60  Title 

Job  Control  Card 

Format  (12,  13,  15,  12,  13,  515,  13,  12,  215,  15.0,  315) 

Columns  3-5  Start  parameter 

= 1 Fresh  set  of  data 

Columiw  6-10  Stop  parameter 

= 1 Stop  after  mesh  generation 

“ 33  Punch  mesh  data  on  cards 
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SAMPLE  -PROBLEM!  NCTCHEO  PLATE 


1 1 

1 

1 1 

400 

1 

1 1 

0.0 

0.0 

1 

2 1 

1.000 

0. 

0 
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4 1 

2.000 

0. 

0 

1 

8 1 

3.000 

0. 

0 

1 

10  1 
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0. 

0 

1 

11  1 
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0 

1 

1 1 

0.0 

0. 

0 

-1 

1 3 

0.0 

2. 
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1 

1 11 

0.0 

4. 
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0.0 
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1 

-2 

11 
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6 

6 7 

1 

3 

7 

11 

12 
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a 
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4 

6 

4 

3 

3 

11 

S 

3 

12 

7 

3 

12 

9 

3 

11 

11 

S 

3 

1 

1 

11 

1 

12 

gure  26 

. Listing  of 

' Input 

to  Grid  Generation  Program  for 

Notched  Plate 
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Columns  21-25  Finite  element  mesh  generation  parameter 

0 Generate  mesh  using  line  segment  data 
= 0 Mesh  data  read  from  cards 

Columns  26-30  Temperature  field 

= 0 Must  be  set  equal  to  zero 

Columns  31-34  Nximber  of  materials  (6  maKimum) 

These  are  the  only  parameters  needed  to  generate  a mesh,  the  rest  can 
be  set  equal  to  zero. 

Mesh  Generation  Card 
Format  (415,  2F10.0,  215) 

Columns  1-5  Number  of  line  segment  cards.  The  actual  number 
of  cards  may  be  less  than  this,  as  long  as  a line 
segment  terminator  card  follows  them. 

Columns  11-15  Number  of  material  assignment  cards 

Line  Segment  Cards 

They  describe  the  boundary  of  the  model. 

Format  (213,  2F8.3,  15) 

Columns  1-3  12  (or  II  if  this  Is  the  first  card) 

1 - index  of  point 

Columns  4-6  J2  (or  J1  if  first  card) 

J - index  of  point 

Columns  7-14  R2  R coordinate  of  point 
Columns  15-22  Z2  Z coordinate  of  point 
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Columns  45-49  IPTION  Integer  value  that  determines  type  of  line 

segment 

Description  of  Line  Segment 

Jump  to  this  point  from  last,  no  mesh  data  will  be 
generated  and  no  line  will  coxmect  the  points. 

Same  as  -1  except  a line  will  be  drawn  between  points. 
Connect  last  point  specified  on  preceding  card  to 
this  one  and  interpolate  into  equal  parts  as 
specified  by  I,  J data. 

Line  Segment  Terminator  Card 
Format  (13) 

Columns  2-3  « -1  Signals  end  of  line  segment  cards.  Must  be 

included  if  number  of  line  segment  cards  is 
less  than  number  specified  on  mesh  genera- 
tion card. 

B -2  The  program  was  modified  to  fix  the  nodal 
points  of  each  vertical  line  in  the  same 
position,  same  Z coord iixate,  as  specified 
on  the  external  boundary  line. 

Internal  Nodal  Point  Generation 

If  -2  is  specified  on  the  line  segment  terminator  card,  the 
program  has  been  modified  to  fix  the  Z coordinate  of  the  internal  nodal 
points  to  those  specified  on  the  external  boundary.  If  this  option  is 


IPTION 


-1 


0 

1 
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specified,  it  is  only  necessary  to  specify  the  positions  of  nodal  points 
on  two  external  boundaries.  For  example,  the  boundaries  with  the  Z and 
R coordinates  of  zero. 


First  Card 
Columns  1-5 

Columns  6-10 

Columns  11-15 

Columns  16-20 

Columns  21-25 
Columns  26-30 


Format  (1615) 

IBMAX-the  I index  of  the  last  column  on  which 
nodal  points  will  be  fixed. 

JBMAX-the  J index  of  the  last  row  on  which 
nodal  points  will  be  fixed. 

NUMJ-the  number  of  J indices  to  be  read  in  on 
the  following  card. 

ICT-I  index  of  the  column  passing  through  the 
center  of  the  crack. 

JCT-JCT  is  set  equal  to  ICT. 

JCMAX-J  index  of  crack  tip. 


Second  Card  Format  (1615)  - This  card  contains  the  J 

indices,  in  sequential  order,  specified  on  the 
boundary  with  Z coordinate  of  zero. 


A Blank  Card  Must  Appear  Here 


Shift  Axis  Card 

If  coordinate  axis  are  to  be  moved  from  the  position  specified  on 
line  segment  card,  the  move  can  be  made  with  this  card.  If  change  is 
not  necessary,  leave  blank.  Format  (2F10.0). 
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Columns  1-10  R shift-shift  in  R coordinate. 

Columns  11-20  Z shift-shift  in  Z coordinate. 

The  following  data  is  needed  to  instruct  the  program  how  to 
perform  the  transition  from  course  to  fine  grid. 

J Omit  Card 

Specified  the  number  of  J rows  which  will  have  segments  omitted. 
Format  (15)  JOMAX 

J Segment  Cards 

Specifies  segment  of  J row  that  is  to  be  included.  Cards  must 
be  in  sequential  order  starting  from  lowest  J-index. 


Format 

(315) 

Columns 

1-5 

JOMIT 

J index  of  'ow  which  segments  are  to 
be  omitted. 

Columns 

6-10 

MSTART 

I index  of  f.rst  nodal  point,  with 
J index  JCMI'K'  that  is  to  be  included 

Columns 

11-15 

MSTOP 

1 index  of  last  nodal  point,  with  J 
index  JOHXT  that  is  to  be  Included. 

I OMIT  Cards 

Specifies  the  number  of  I columns  which  will  have  segments 
omitted. 

Format  (15)  lOHAX 
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I Segment  Cards 

Specifies  segment  of  the  I column  which  is  to  be  included  in  the 
grid  pattern.  Cards  must  be  in  sequential  order  starting  with  lowest 
I index. 


Format  (315) 
Columns  1-5 

Columns  6-10 


Columns  11-15 


lOMIT  I index  of  column  which  segments  are 
to  be  omitted. 

LSTART  J index  of  first  nodal  point,  with 
I index  of  lOMIT  that  is  to  be 
included. 

LSTOP  J index  of  last  nodal  point,  with  I 

index  of  lOMIT  that  is  to  be  included. 


Expansion  Parameter  Card 

The  information  on  this  card  signifies  the  beginning  and  the  end 
of  the  fine  gridding  in  the  crack  tip  region,  in  terms  of  the  J index. 

Format  (315) 

Columns  1-5  lEXPAN  The  J index  of  the  row  following  the 

last  JOMIT  index  specified  in  the  J 
segment  data. 

Columns  6-10  IIEXD  I index  of  a node  point  which  tells 

the  program  that  after  this  point  is 
passed  it  is  to  begin  generating 
triangular  elements  when  the  nodal 
points  are  in  the  correct  formation. 
This  point  is  directly  below  the 
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crack  tip  and  in  the  row  with  a J- 
index  one  less  than  the  first  JOMIT 
index  specified  in  the  J segment  data. 

Coliflnns  11-15  JJEXP  J index  of  the  nodal  point  described 

in  IIEXP. 


Material  Block  Assignment 

Each  card  assigns  a material  number  to  a block  of  elements 
defined  by  I-J  data.  The  number  of  cards  must  agree  with  the  number 
specified  in  columns  11-15  of  mesh  control  data  card. 


Format 

Columns 

(515) 

1-5 

Material  definition  number 

Columns 

6-10 

Minimum  1 

Columns 

11-15 

Maximum  1 

Columns 

16-20 

Minimum  J 

Columns 

21-25 

Maximum  J 
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C SOLIDS  WITH  different  jRTHOTROPIC.  TEMPERATURE-DEPENDENT  MATERIAL 
C PROPERTIES  IN  TENSION  AND  COMPRESSION  INCLUDING  THE  EFFECTS  OF 

C internal  pore  fluid  PRESSUHES  AND  thermal  STRESSES. 

««««»*««»«•***«•«•««*«««««*#«*««*«««***««*»«««««*«»«• 
COMMON/BASIC/NOMNP.NUMEL.NUMPC  NUMSC.ACEL2»ANGVEL«TREF.V0L.IFREQ 
COMM0n/NPDATA/R(100O) .code (1000) «XR (1000) .2(1000 }»XZ (1000). TdOOO) 
DUdBLE  precision  CRZ.XI.1RR.22.S.RRR.222 
COMMON/ELOATA/IX (1000.5) .EPR(lOOO) . ALPHA ( 1000 ) .PST ( 1000 ) 
CQMMON/SOlVE/NCODE  (30.100)  .N'UMTC 

COMMUN/TO/IMIN(lOO) .IMAX(lUO) «0MIN(30) .UMAX (30) .MAXI .MAXU .NMTL. NBC 
1 .MINI.MINJ 
COMMON  /OHG4/  ISTOP 
C0KM0N/LC2/  KP.RP (250) .ZP (250) 

COMMON/ JEH/NCHC.IMIV  (30 ) .IMAV  (30)  .JMIV(30)  .JMAVOO) 

C0MM0N/0HG5/  JS.TOP 
INTEGER  TITLE  (20) 

c read  and  write  inplt  Data  card  images  for  all  cases 

OSTOPsU 
200  CONTINUE 
CALL  DATA 
ZFORCt  = 0,0 
ELMASS  » 0,0 
VOLM-0.0 

C ««•««•*« ««•*»»•**«  •«««•«« •«••«•«««• 

C READ  AND  WRITE  CONTROL  INFORMATION 
C 

READ («» 1000)  title. NPP.ISTART. ISTOP, lOEF.IPLOT.NNLA.lMESH.NUMTC. 

1 NPOHPR.iNUMNAT,NUMPO,NUMSC.TReF,NTCA,lFHEO 

2 .IBPLT.NUMNP.N'JMEL 

WRITE (0, 2000) (TITLE (I). 1=1. IS). 

1 ISTART. ISTOP. IDCF.IPlOT.NNLA.IMESH. 

1 NUHTC.NPORPR.iNUMHAT.NUMPt.NUMSC.TBEF 

2 .lOPLT 

C OBNERATE  FINITE  ELEMENT  MESH 

215  CALL  MESH 

C 

C INITIALIZE 

C 

C HEAD  ANO  WHITE  NOOAL  POINT  AND  ELEMENT  DATA 

2£5  call  points 

C •••••••• w*************************************^************************ 

C OUTPUT  ELEMBNT  OATa 

PRINT  3100 
HPR1NT"0 

00  350  Nwl.NLHEL 
|F(MPHINT,NE*0)GO  TO  300 
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WR1TE(6*200<3) 

MPRINTbSO 

300  MRBINT*MPRlNT-i 
NSTRSS>20 
THICKsO.281 

350  >«RITE{beg009)  Ni (IX (N» I J , I»1 15> #ALPHA <N>  »T(N) ,PST(N) .NSTRSStTHICK 
IF(  ISTOP  .N£,  33  ) 60  TO  362 

00  360  NPUN>1«NUMEL 

PUNCH  3002«NPUN( (IX(NPUN*I) tlaltS) tNSTRSStTHICK 

360  CONTINUE 

00  361  NP>l«iNUHNP 

PUNCH  3200tNP*R(NP) tZiNP) 

361  CONTINUE 

362  CONTINUE 

370  IFdSTUP.EOtl. OR. ISTOP. EQ. 11)60  TO  910 
910  1F(  NPP  .NE.  0 ) 6C  to  200 

IF  ( 16T0P  .EQ.  1 ) 60  TO  200 
IF  ( IFREQ  .EQ.  0 ) 60  TO  200 

C SUN  THE  ELEMENTS  OF  THE  Z-OIRECHON  BOUNDARY  FORCE  MATRIX 
00  905  IFOR  ■ItNUMNP 
90S  ZFCRCE  a ZFORCE  « XZ(IFOR) 

ZFOHCE  « ZFORCE  • 2.  • 3.IA15927 
ELMASS  m EL'FASS/2.0 
WEIGHT  « ELFASS*12.0«32.17 
QFORCt  ■ EUFASS«ACELZ 
ELMASS  • EUFASS*I2.0 
WRITE  (6*3000)  ELMASS 
WRITE  (6*3001)  ZFORCE 
WRITE(6*3003)WE1QPT 
WR1T£(6*3004)6FORC£ 

60  TO  200 
C 

1000  FORMAT(20A«*/» 

1 12*I3«t5«I2*I3*315*I3*I2*2l5*F5. 0*615) 

1001  format  (AFIO.O) 

2000  F0RMAT(IH1»1SAA*/* 

1 33H0  START  PARAMETER— *14  */ 

2 33HQ  STOP  FARA“tTER“«»»"““-*'*’““»’*“'*-* I* */ 

3 33M0  IF  1,  PLOT  DEFLECTIONS  *I4»/ 

A 32H0  IF  1*  SMALL  PLOT,  IF  g.  LARGE*/ 

5 33M0  PLOT.  OTHERWISE  NO  PLOT.— — *-» I**/ 

6 33H0  NUMBER  OF  APPROXIMATIONS-— 14*/ 

7 33M0  IF  1*  GENERATE  MESH— — — — »I4*/ 
a 33H0  NUMBER  OF  TEMPERATURE  CAROS— *14./ 

3 33H0  NUMBER  OF  MATERIALS— ——— *14*/ 

A 33H0  NUMBER  OF  EXTERNAL  PRESSURES— *14*/ 

5 33HQ  number  of  SHEAR  CAROS— ——*»  14*/ 

6 33M0  REFERENCE  TEMPERATURE— —«»E12. 4./ 

7 33H0  boundary  PUT  OPTION  .14*//) 

2001  FORMAT  (7TH  A FUNOAMENTAL  FREQUENCY  WILL  BE  COMPUTED,  A LONGER  RUN 

1 TIME  WILL  0£  OBSERVEO/BOH  OUE  TO  THE  NEED  TO  RECOMPUTE  EACH  ELEME 
2NT  STIFFNESS  MATRIX  IN  SUBHOUTINE  STRESS) 

2002  FORMAT  (24H  THE  ANGULAR  VELOCITY  1S,E12.4*/31N  AND  THE  AXIAL  ACCEL 
lEHATlON  IS  .EI2.4) 

2003  FORMAT  (2SN  THE  R ACCELERATION  Is  »E12.4/27H  ANO  THE  Z ACCELERATIO 
IN  IS  *LI2.a) 


71 


2004  FORMAT  (//42h  THE  PLANE  STRAIN  OPTION  HAS  BEEN  SELECTED) 

2005  FORMAT  (//4ap  THE  PLANE  STHESS  OPTION  HAS  BEEN  SELECTED) 

2006  FORMAT  (BeHl  EL  I J K L MATERIAL  ANGLE  TEMPERATURE 
I PRESSURE  PRINT  THICKNESS  ) 

2009  F0RMAT(I5,4l4f I8.Fll.l.2F13.3*5X,I5,5X»F10,4) 

2013  FORMAT  (30Hl  PRESSURE  BOUNDARY  CONDITIONS) 

2015  FORMAT  (27HI  SHEAR  BOUNDARY  CONDITIONS) 

2020  FORMAT  (/29H  THE  PROCEDURE  CONVERGED  IN  «I2*34H  TENSION  . COMPRES 
ISICN  ITERATKNS  ) 

2021  FORMAT  4/36H  THE  PROCEDURE  DID  NOT  CONVERGE  IN  tI2.33M  TENSION  - 
ICQPPREbSlON  ITERATIONS) 

2022  FORMAT  (/29H  THE  PROCEDURE  CONVERGED  IN  fiataOH  NONLINEAR  ELASTIC 
I ITERATIONS  ) 

2023  FORMAT  (/36H  THE  PROCEDURE  DID  NOT  CONVERGE  IN  ,I2»30H  NONLINEAR 
lELASTIC  ITERATIONS  ) 

2030  format  (/51H  NUM8ER  OF  lENSION-COMPRESSION  APPROXIMATIONS----* 

I 14//) 

3000  FORMATdM  «//*10X*  4HMASS.20  (2H.  , ) *E20 . 14 . 3X.5HSLU6S) 

:001  FORMATdH  »/»lOX.3SMAXIAL  FORCE  DUE  TO  NORMAL  STRESSES./. lOX. 

1 40HSHEAR  stresses.  AND  CONCENTRATED  RADIAL.  /. 

21H  .lOX.  UHAND  AXIAL  FORCES  *14(2H.,)*  E20.H.3X.3HLBF) 

3003  FQHMATdH  ./.lOX.  6HWEIGHT.  19  (2M.  , ) , E20.14.3X.  3HL8F) 

3004  FORMATdH  ./.10X.20HAX1AL  INERTIAL  FORCE  . 12  (2H.  . > .E20 . 14.3X. 

I 3HLBF) 

3002  FOHMAT(6I5.2CX.I5.SX.F10.4) 

3100  FORMATdH  . 12MELEMENT  DATA) 

3200  F0RMATd5.40X.2F10, 4) 

ENC- 
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C MESH—  THIS  ROUTINE  GENERATES  THE  MESH  GIVEN  LINE  SEGMENT  INPUT 

subroutine  mesh 

c««***«*»«»«««*««*«*««*«««**«««««****«««*««*«*«*«*#««*«««*«»*«***«*«*** 

C NCOCE(I.J)  = 1 FOR  BOUNDARY  DEFINITION  POINTS 
C > 2 FOR  INTERPOLATED  BOUNDARY  POINTS 

C *2  FOR  INTERPOLATED  INTERIOR  POINTS 

C * A FOR  EXTERIOR  POINTS 

C a .£  FOR  VOID  POINTS 

C IPTION  a 0 FCR  SINGLE  POINTS 
C * I F'CR  STRAIGHT  LINES 

C B 2 F'CR  INTERNAL  DIAGONAL 

C *3  FCR  3-PCINT  ARC 

e a A FCR  2 POINT  ♦ CENTER  ARC 

C a 5 pcR  2 POINT  ♦ RADIUS  ARC { INITIALIZATION  ONLY) 

C a 6 FCR  2 POINT  * RADIUS  ARC 

(;»«««•«*•»««•«•««««••«•*«*•«•««»«»•*•*«««*••««««•**•«••«««•#*•«•««»«««* 

C 

COMMON/TO/IMIN(100)  tIHAXaOO)  tUMINOO)  •JMAXOO)  ^MAXI  »MAXU*NMTL*NBC 
1 •MINItMl'NJ 

COMMON/NPOATA/R  (1000)  *COOE(  1000)  tXR(  1000)  .ZaOCO)  1X4  U 000)  tTdOOO) 
CQMMON/€LDATA/IXaOOOtS)  .EPR  (1000)  ♦ ALPHA  (1000)  »PST  { 1000) 
COMMON/UC2/  KPiRP(250) .ZP(250) 

COMMUN/6ASIC/NUMNPtNUMEL»NUMPC.NUMSC*ACELZ»ANev£L»TR£F.VOLf IFREO 
COMMON/OHGA/IaTOP 

CUMMON/U£H/NCHC«IMIV(30)  f IMAV(30)  •jMIV(30)  tJMAVOO) 
COMMON/SOLVE/NCOOE (30.100) .NUMTC 
DIMENSION  AR(30.100).A2(30.100) 

DIMENSION  JUG (20) 

equivalence  (ru)  .aru.d  ) t (Zii)  lAza.i)) 

C MESH  CUNTHOL  INFORMATIUN 

HEAD  (V»iOOO)NSEG.*N8C»NKTL»NLlM.CONI»CONo.iSET»USET 
mRITE  (6.2000)  NS£6.<KBC. NHTL.NLIM. CONI  tCQNu.lSET.vJSET 

C initialize 

Kpal 

ISEOaO 

PI«3.lAXS927 

n»i 

IC«0 

DO  110  U«1.100 
00  iOO  I«1.30 
NCC0t(l.J)a4 
AR(I.O)aO, 

100  AZ(l»0)a0. 

lMlN(J)a3o 
UU  IMAX(U)aO 

DO  120  I«1.30 
UHAX(I)aO 
120  UHIN(l)alOO 

C LINE  SLGMENT  CAROS 

PRINT  2001 


NEAD(9«1001)12tJ2*P2?Z2 

IPTION*! 

13U  PRINT  2005 

PRINT  2010fI2t^2tR2*22 
200  ISEQ*Ii»E84>I 
11*  IS 
01*02 
Rl*R2 
21*22 

AR(U«01)»R1 

AZ(I1,01)*21 

RP(KP)»R1 

IP(IPTION,LT.O)RP(ttP>*-Hl 

2P(KP)«21 

KP*KP«1 

NOUOE(11»OU«1 
CALL  HNlHXdltOl) 
lF(ISE8.Le.N§EQ)  6C  TO  500 
IFtlC  .EQ.  I ) 60  TO  2A9 
250  RtAO{9tl00l! I2t02«R2*22f I3t03«R3fZ3«lPY10N 
249  IF(I2«NE.  >2«ANO«lC.EQ.O)  60  TO  251 
IF(I2.NE.  -2)  60  TG  252 

READ(9>5000i  IBRAXi JBHAX*NUMO«ICT*OCT*oCMAX 
SOOU  FORMAT (1615) 

HEAD(9tSOOO)  (OOBOX) «N*1«NUM0) 
ii»n*i 
00«1 
i2«n 

02»00B(U 

H2*AR(12«1) 

22«A2a,  02) 

IPTION*”.! 

I3»0 

03*0 

R3*0.0 

23*0.0 

IC*1 

60  TO  OSi 

252  IF (12  «6T.  18MAX  ) GO  TO  255 
IF(02.bQ.0BMAX)  00  TO  253 
00*00*1 

12*1! 

02*008(00) 

IFa2.0E.ICT.AN0«12*LE.OCT»ANO.O2.GT.OCHAX)  60  TO  253 
IPTION*! 

60  TO  254 

253  n*ii*i 
00*1 

i2*n 

IF(I2.6T*  IQPAX)  OC  TO  255 
02*008(1) 

IPT10N*»1 

254  R2*AH(12«1) 

22«AZ(li02) 

IC*1 

60  TO  251 


o o o o o o 
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255  IC*0 
I2a-i 
60  TO  251 

251  IF(I2  .EO.  -1  ) GO  TO  500 
IF(IPTION,LT.O)GO  TO  130 

WRITE (b» 20 10) I2»J2tR2iZ2«I3» J3.R3»Z3,IPTI0N 
IPTION=IPTION-H 

60  TO  ( 200«300»300.346t3^6,346.346)  ♦ IPTION 

C ««»«««««  ]«>«*«««««««««««««««««*«««*«««««««««« 

C 6BNERATE  STRAI6HT  LINES  ON  BOUNDARY 

C***«**4'  ««»««**««»«««««««««*««•»»««<»««*»«««*«*««•««««««««««««»«««««««««« 

300  0l8AbS(FL0AT  (12-11)  ) 

OJ»ABS(FLOAT(02-J1) ) 

ISTRTsIl 
ISTP8I2 
USTKTsUl 
JSTP=02 

UIFFsAMAXl (OlfOd) 

ITERaOIFF-l. 

IINC»0 

dl'NCsO 

IF'I2.NE.I1)  IINC«(I2-I1)/IABS(I2-I1) 
lF(J2.NE.sJl)  dINCa(J2-Ul)/lABS(U2-dl) 

KAPPAsl 

lF(I2,NE.n.ANO.J2.NE.dl.AND.IPT10N.NE.3)  KAPPA»2 
IF(KAPPA.E(J.2)  0IFF»2.*0IFF 
RINC«(H2-R1)/0IFF 
Z1NC=(Z2-Z1)/0IFF 

WRITE  (b, 2002)  Dlv<DwtDlFFtHINC»2INC.ITER#IlNC»JINC*KAPPA 

CHECK  FOR  input  ERROR 

IF(IPTION,E(J.3  .ANO.OI.NE.OJ)  60  TO  310 
IF (KAPHA.NE.e.OR.OI.EU.DJ)  60  TO  320 
310  WRITE(bt2003) 

60  TO  200 

INTERPOLATE 

320  I>I1 

1 

WRITE{bf.2004) 

WRITE(b»2004) 

00  340  H»1«ITER 

IFdTEH.EO.O. AND. IPTION. £0.2)60  TO  345 
IF(KAPPA.EOr<)  60  TO  330 
lOLOsI 
I»1+1INC 
JOuOaJ 
JWw^'UINC 

AR(If  J)»AR(KLD*vJQLD)*HINC 
AZ(IfU)«AZ(I'CLO*dOLO)*ZINC 
WRITE  (6»2005)  I*JfAR(I»0)  »AZa«J) 

CALL  MNIMX(I*J) 

IF(NCOOE(Ifd) .NE,4iWRlTE(6»2100) I«d 
N0C0E(I«d)32 
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60  TO  340 
330  IOUt)*I 
I*I+IINC 

AR(ItJ)3AR(ICL0eJ)4RINC 
AZ(ItJ)3AZ(KL0tJ)*ZlNC 
MRn&(bi2005)  ItUiARCItJ)  lAZdt  J) 

IF (NC0U£(I«3) .NE.4)WRlTE(6i2100) If J 
NOCDEdf  J)s2 
CALL  MNIMXdfJ) 

OQLOej 

Jay+JINC 

ARd»J)«AR(I.J0LDJ4RINC 
AZdf  J)sA2(IiJOLO)*ZINC 
IF  (NCOUEdfJ)  .NE.4)  WRITE  (6f2I00)  IfU 
NO0DE(IfO)32 

WRITE (6.2005)  I . J. AR (I . J) . AZ (I » J) 

CALI  MNIMXd.J) 

340  COldINUE 

345  IF (KAPHA. EG. i)  60  TO  200 
IULDkI 
I*I*IINC 

AR(I.J)3AHd'CLU.J)fRINC 
A2d»J)3AZ(KL0f  J)*ZINC 
IF(NCOUEd.J)  .NE,4)WRITE(6»2I00)  Ifd 

NOCOE(I»J)32 

WRITE  (6.2005)  I.JfARdtJ)  .AZd.J) 

CALL  MNIMXdiJ) 

60  TO  200 

346  WRITE(6fA00i) 

4001  FORMAT!  *0«  »•**«*  INPUT  ERROR  —.VALUE  OF  IPTION  GREATER  THEN  I 
I— — -EAECUTKN  TERMNATEDd 
STCP 

500  MAAI3U 
MAAOsU 
NIiM«30 
MIINU3I00 
00  503  I«1.30 

IF(UMAXd)  ,6T.MAXU)MAXvi3jMAXd) 

IF(UMlNd)  ,LT.MINvi)MINU3jMINd) 

503  CONTINUE 

00  507  UsltlOO 

IFdMAX(U)  .0T.MAXI)HAXl3lMAX(J) 

IFdMIN(U)  .LT*MINI)MINI«IHIN(U) 

507  CONTINUE 

READ  (VflOOO)NOHC 
IF(NOHC.E(9.0)GO  TO  511 
HRlTE(6f6500) 

00  508  NOsliiNOHC 

READ  (V.8000) INI. IMAfUMIf JMA 

WRITE (6. 8000 )IMIf£MA*UMI*JHA 

1HIV(NU)3IMI 

lMAV(NO)aIMA 

JM1V(N0)3JNI 

UMAV(NU)«OMA 

00  508  JsjMIiJMA 
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DC  508 

508  N0CDE(I»J)s5 
8000  F0KMAT(415) 

8500  F08MATUHO.  ZOhVOlD  SPECIFICATIONS  ,/t 
1 «lh  IFI  IMA  JMI  JMA  ) 

IF(NUHC,LE.30)GO  TC  511 
WRITE (8*9000) 

9000  format (33HO***«ERHCR****TOO  MANY  VOID  CAROS) 
STCP 


C CALCULATE  COORDINATES  OF  INTERIOR  POINTS 
(;»***»■»********«#*«  ***#***»***********, 

511  IF  ( MAXI*1  ,LE.  2 ) 60  TO  530 
IlaMINUl 
I2»MAXI-1 

IF(NLIM.LT.l)  MIM  s IQO 
DO  520  N=1,NLIM 
HESIDsU. 

DO  510  1*11*12 
JlaJMINd)*! 

U2»JMAA(I)-1 
DO  510  Uaol*.w2 
KOCE»NCODE (I*U) 

60  T0{bl0«610t509»506«510) *KODE 
506  NOCDE(I.J)a3 

509  DHs(AR(I*1,J)*AR(I-1,U)*AR(I,J+1)*AR(I.J-1) )/4,-AR(I.U) 

1 *CONI  * (AR(I*1,U)  - AR(I-1,U))/FL0AT(8*(I*ISET) ) 

2 ♦ CONO  * (Af>(I*J»ly  - AR(I*J-1))  / FL0AT(8*(J*USET)  ) 
DZa<AZ(I*l,J)*AZ<I-l*U)+AZ(I*J*l)*AZ(I*J-l) )/4.-AZ(I*J) 

1 ♦ CONU  * (AZdiJ^l)  • AZd»U-l)  )/FL0AT  (8*(U*USET)  ) 

2 * CONI  * (AZd*lfw)  - AZd-l*U))  / FL0AT(8*d*lSET)  ) 
RBSIDsRESIO^ABS (DR) +ABS(DZ) 

AR(I. J)»AR (If J) ♦1,8«DK 
AZd.J)»AZd*J)*l,8*DZ 

510  CONTINUE 
IF(N,EU,1)  RES1»RESID 

IF  (N.EU.I.ANC.RESIC.EU.O)  60  TO  530 
IF(RESID/RES,\,LT.1*E-4)  60  TO  530 
520  CONTINUE 
530  WRITE (6*2009)  N 
KPaKP-1 
600  CONTINUE 
999  WRITE(6»4000) 

4000  FORMATdH  * 9HEN0  MESH  ) 

C #«**##**  *)###**#*(*****>****  ************** 

0 ****##***•###**###*####******#******#******************************  ******************************^****^^^^**** 

c 


1000  F0HMAT(4I5*  2F10, 0*315) 

1001  format  (2(213*2F8.3)*I5) 

2000  FORMAT  (30HI  MESH  GENERATION  INFORMATION// 

3 41H0  NUMBER  OF  LINE  SEGMENT  CARDS———— 

4 41H0  NUMBER  OF  BCUNDARY  CONDITION  CARDS— 

5 41H0  NUMBER  OF  MATERIAL  BLOCK  CARDS—— 

4 41H0  NUMBER  OF  BCUNDARY  CONDITION  CARO«— 

5 41H0  NUMBER  OF  MATERIAL  BLOCK  CAROS—— 

6 41H0  NUMBER  OF  ITERATIONS— 


.13*/ 

*13./ 

*13*/ 

*13*/ 

*13*/ 

*13*/ 
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7 41H0  POLAR  coordinate  PARAMETER  I=»— — — *E12.4,/ 
tt  41H0  POLAR  COORDINATE  PARAMETER  J— — — *E12.4*/ 

9 41H0  I CURVATURE  MODIFICATION— I3t/ 

I 41MO  J CURVATURE  MODIFICATION— 13*///) 

2001  FOhMATdHO.  eih  INPUT  I J R Z 12  J2  R2 

1 Z IPTION  > 


2002  FORMATUh  , 5H  DIa*F4.0*  5H  DU»*F4.0t  7H  DIFF»*F4.0. 

1 7h  ring*. Fa. 3.  7H  ZINC»*F8,3*7H  ITER«*I3. 

2 7H  IINC».I3.  7H  JlNC*»I.i,  8H  KAPPA*. 12) 

2003  FORMAT  (1)(.38H**BAD  INPUT— T::aS  LINE  IS  NOT  DIAGONAL) 

2004  format  <30R  I J AR  AZ) 

2005  FORMAT  (2I5.2F10.3) 

2006  format  (51H  «*  BAD  INPUT  - THESE  POINTS  DO  NOT  DEFINE  A CIRCLE*/* 
13X*<bFU.A,10)t*2E20,a) 

2007  FORMATUH  , 2IH  CENTER  COORDINATE  { *F8 ,3  . IH*  *F8.3  . IH)  ) 

2008  FORMATUH  , 7H  ANQl«,F9.6*7H  AN62**F9.6*7H  DIFF«*F3,0* 

1 9H  DELPHI*. F9. 6) 

2009  FOHMATUhO.  30H  CCORDiNATtS  CALCULATED  AFTER  *13* 

I IlH  ITERATIONS  ) 

2010  F0RMAT(7X, 2(214, 2F8.3)*I6) 

2100  FORMAT  (54h****«WARMNG******  NODAL  POINT  KITH  (I,J)  COORDINATES  (* 
1 I2*1H,,I2,21H)  HAS  BEEN  RE-DEFINED) 

3000  F0hMATU6I5) 

3100  F0HMAT(2El5.e,II0) 

3200  FORMAT(8F10*'5) 

RETURN 

END 
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C 

C 


POINTS—  THIS  HCUTINE  ASSIGNS  MATERIALS*  TEMPEHATUHESt  ETC. 

subroutine  points 


I 


COPMON/BASIC/NUMNP«NUMEL*NUMPC»NUMSC*ACELZ*AN6VEL»TREP»VOL» IPREQ 
COPMON/NPDAT.^/R  (1000)  tCOOE  { 1000 > *XR  < 1000 ) »Z  < 1000)  »XZ  { 1000)  tT  ( 1000) 
COPMON/ELOATA/IX (1000*5) *EPR(1000) *ALPHA(1000) *PST(1000) 

DOUBLE  PRECISION  X*Y*TEM 
CQNMON/SOLVE/NCODE (30*100) *NUMTC 


COPMON/TD/IMIN(lOO) »IMAX(100) *JMIN(30) *JMAX(30) *MAXI *MAXU,NMTi  *NBC 

*mini,minj 

COPMON/PLANE/NPP 


COKMON/JEH/NCHC*IMIV(30) *IMAV(30) *UM1V(30) *UMAV(30) 
C0PM0N/0HS4/  ISTOP 

DIMENSION  AR(30*100) *AZ(30»100) *MATRIL (20*5) *BLKAN6(20) 
OIMENSION  IOROER(30*100) 

DIMENSION  JOMIT(AO) *MSTART(40) ,HSTOP(AO) 

DIMENSION  lOMIT(AO) *LSTAHT(40) *LSTOP(40) 

DATA  IUROER/3000«0/ 

EQUIVALENCE  <R ( 1 ) * AR ( 1 , 1 ) ) * (2 ( 1 ) . A2 ( 1 » I > ) 

C ***«#***##**#*» 

READ(9.4999)  RSHIF*ZSHIF 

4999  FOHMAT(2F10.0)  ^ 

HEAO(9.5000)  UOMAX 

5000  F0RMAT(I5) 

00  7 UO«l,dOMAX 

READ (9*5001 ) dOMlT (JO) *MSTART (dO) *MSTOP (dO) 

5001  F0RMAT(3I5) 

7 CONTINUE 

H£AO(9*5004)  IOMAX 
5004  FOHMATdS) 

DO  6 lUal, IOMAX 

6 READ(9*500I)  lOMIT  (10) *LSTART (10) *LSTOP (10) 

HEAO(9*5001)  I£XPAN*IlEXP*ddEXP 
NPsO 


dOai 

DO  120  daMINwtMAXd 

IF(  d .EO.  dCMIT(dO)  ) GO  TO  90 

NSlARTslMlNOv) 

NSTOPalMAX (d) 

60  TO  91 

90  NSTART-MSTART  (00) 
NSTOPbMSTOPOvO) 
d0sd0*l 

91  10«1 

DO  120  IbNSTART*NS70P 
IOROER(I«d)«0 

IF(NC0UE(I*d),E(J,4)  GO  TO  120 
IFd.NE.IOMlTdO))  GO  TO  96 
IF(d.Lt.LSTOPdO)  ) GO  TO  92 
60  TO  93 

92  lF(d  .GE.  LSTARTdO)  ) GO  TO  95 

93  I0aX0«I 
60  TO  120 

95  10«I0*1 

96  NPaNP^l 
RCNP)BANd*d) 


Z(.NP)sAZ(I.J) 

10RDER(I»J)«NP 
lEO  CONTINUE 
NUPNPsNP 

C READ  AND  ASSIGN  BOUNDARY  C0.4DITIONS 

£«««««*•«*»«*«»*«**«*•«*««**«<»««**•«««*««*«««**«****»**«•«*«««««««»*««**« 

C INITIALIZE 

C 

DO  150  iBitlOOO 
T(I)=0.0 
ISO  PST(I)eO.O 

DO  ZOO  I»1«NUMNP 
COD£(I)sO. 

IP(R(I) .EQ.0..AND»NPP.EQ*0$  COOE(I)  ■ 1. 

XR(I)sO. 

XZ(I)rO. 

200  CONTINUE 
C 

IF{NUC«EU.0)  GO  TO  220 
00  210  IbCONilfNBC 

READ (9 t 1002)  II»I2«Jlf J2«CUN«RC0NtZC0N 
00  210  Iall,I2 
DO  210  JaulfwZ 
NP«IOROER(I*is.) 

COCE(NH)aCON 
XR(NP)sHCON 
210  XZ;NP)sZCON 
220  MPHINTaO 

PRINT  1300#NIMNP 

NP«0 

JO«I 

DO  240  JaHINv«NAXJ 

IF(  U .EO.  OOMIT(v)O)  ) GO  TO  225 

NSTARTalMINPw) 

NSTUP«IMAX(J) 

60  TO  226 

225  NSTARTaMSTART(OO) 

NST0P»MST0P«-.0) 

UOa>JO«l 

226  10«I 

00  240  I«NSTART.NS70P 
IF  ( NCOOEdtJ)  .EG.  4 ) GO  TO  240 
IFd.Nb.lOHIT(IO)  ) GO  TO  256 
IF(J  .LE.USTCPdO)  ) GO  TO  254 
GO  TO  257 

254  IF(U  .GE.  LSTARTdO)  ) GO  TO  255 

257  I0»I0*1 

GO  TO  240 

255  10*10*1 

256  NP*NP*1 

IF(MPRINT.NE.O)  GO  TO  230 
NR1TE(6»2000) 

MPRINT*50 

230  MRRINT*MPRINT-1 
R(iNP)*H(NP)«fiSNIF 
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2(iNP)xZ{NP)-2SMIF 

WRlTt(b»200n ltJtNPtC00E(NP)«R(NPi«Z(NP) «XR(NP) «XZ(NP) «NC0DE(ItJ) 
240  CONTINUE 

C ASSION  MATERIALS  LN  BLOCKS 

00  300  MlaltlOOO 
300  UJMl,a)»0 

00  310  IMTLBitNHTL 

REAO  (^*1000)  MTL» (MA7RIL(IMTLtIM),IM«2.5)»BLKAN6(MTL) 

310  MATRILdHTLtDxMTL 

C ESTABLISH  ELEMENT  INFORMATION 

«««««**« 

UJMAX«MAXo-l 

lOCOEal 

JOCDEal 

JvJCOOEal 

N£L>0 

JUaO 

IlaO 

60  TO  312 

311  CONTINUE 

312  JJaJU«l 

1F(JJ  .6T«  UvMAX)  .CO  TO  400 
UUC00E«1 
11*0 

NSTANT«IMINOvO> 

NSTOPalMAX(Jw) 

IlaNSTART 
QC  TO  31h 

313  II«II*1 

iFdl.OE.NSKP)  60  TO  311 

314  IF(IUROERai*jU)  ,EQ.  0 ) 60  TO  313 
IIP»II+1 

UUP*vJ  1 

IF  ( NCUQEdltUU)  .EO.  4 ) 60  TO  317 
lF(NCOUEdlPtJU  ).EQ.A)60  TO  317 
IFlNCOUEdl  »Js)P»  .£0.4)60  TO  317 
XFlNCOOEdlPiJOP)  .£0.4)60  TO  317 
1F(JJ.6£.OUEXP.ANO.II«6E.IIEXP)  00  TO  316 
IFtNCOOEdI  »JU  )»N£.5)60  TO  316 
IF(NCOOEdlP.JvJ  )UN£.5)6C  TO  318 
lF(NCODEdI  iJuP)  .iNE.b)80  TO  316 
lF(NCOUEdlPf OOP)  .iNE.5)60  TO  318 
DO  315  NOalflNOHC 

lFdI.6E.IMIV  (NO)  .AND.II.LE.IMAV(N0).AND. 

1 UU.6E.UMIV (NO) .AN0.J0.LE«0MAV(N0),AN0. 

2 1IP.6E,IMIV(N0) .aNO.IIP.LE.IMAV(NO) . AND. JJP.OE .UMIV (NO) .AND. 

3 0UP.LE.uMAV(NO))60  TO  313 

315  CONTINUE 

316  JvJCOOE>2 
60  TO  318 

317  UvJC0UE«2 
60  TO  313 

318  CONTINUE 
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C CHECK  FOR  DIHECTICN  OF  EXPANSION 
IF(JJCUDE,EQ.2)  60  TO  320 

319  CONTINUE 

IF(I0HUER(II*JUJ.GT,0. AND. lORDERdl. OOP). GT. 0)60  TO  333 
IF(IORLlER(IIiOU).GT.O.AND.IORDER(lIPfJj).GT.O)  60  TO  321 
WRITE(6i6000)  II«Uu»IIPtUJP 

bOOO  FORMAT)*  • » * **n*ERRCR****»  POINTS  * *4I5»  * DO  NOT  DEFINE  A VALID 
* EXPANSION  EXECUTION  TERMINATED*) 

LTQP 

320  CONTINUE 

IF(  JJ  .GE.  IEXPaN  ) 60  TO  319 
1F(I0RUER(IIP<0J)  .EQ.  0 ) 60  TO  313 
lF(IORDER(irP«UJP).EQ.O)  60  TO  321 
60  TO  319 

321  continue 

C EXPANSION  IN  U’DIRECTION 

IFtIORUER(IIfUuP).EQ.O.AND.IORDER(lIP»UUP) .EQ.O)  60  TO  331 

IF (lORUER(II.OUP) ,6T.O,AND.IORDER(IIP»dJP) .6T.0)  GO  TO  329 

IF(I0RDER(II«JUP) .EQ.O)  60  TO  322 

IITEMPall 

JdTEMPsJdP 

JCCDEnd 

60  TO  331 

322  UTEMPellP 
OUTEMPojJp 

J000Ea2 
60  TO  331 

323  CONTINUE 

c triangular  elements 

lAsIORUERdItJU) 

IBalOHUERdIPfUU) 

ICsIOROERdITEMP,dvTEMP) 

lOalC 

60  TO  332 

324  lAalORUERdliJU) 

IB«10R0ERdITEMP,dgTEMP) 

IC«10R0ERdl«vJUP> 

lOwJC 
U0C0Eb4 
60  TO  332 

325  lAsIOROERdItJUP) 
l45BlORUERdITEMP,J„TEMP) 

ICalOROERdlPfOOP) 

lOalC 

JCCDEsl 
GO  TO  332 

326  lA«IOROERdI*JU) 

IBalO.ROERdlPiUJ) 

ICalOROF-RdlTEMPtOuTEMP) 

IO*IC 

60  TO  332 

327  1A»I0R0ER{HF»UU) 

IBalOROERdlPtUUP) 

ICalOHDEHdneMPfdwTEMP) 

IDalC 

dOODEab 
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I 


j 


60  TO  332 

326  IA3lOHOER(ItP«OJP) 

16>I0RUER(IIi JUP) 

IC«IORUER(HTEMP»JwTEMP) 

IO»IC 
JOCUEsl 
60  TO  332 

329  60  TO  (330t323«326)  tJCOOE 

330  1A=I0R0ER(II» 03} 
lB>IORDER(IIPtOO) 

ICsIORDERdIPiOOP) 

IDsIORUER(IliOuP) 

60  TO  332 

331  OOP*OjP«1 

IP(JJP.GT.MA>0)  CALL  ERROR U lOOPtMAXOi II ) 

60  TO  321 

332  NEl.BNEL*! 

IP(NEL.6T.1000)  call  ERRORv2.NELi1000«0) 

1X(NEL*1)»IA 

IX(NEL«2)>IB 

IX(NEL»3I*IC 

IX(N£Lt4)«I0 

IX(NEL>5)>MTL 

ALPHA  (NED  oBLKANG  OPTD 

60  TO  ( 345»324«327 1325*326  ) tOCOOE 

333  CONTINUE 

C EXPANSION  IN  I-OIRECTION 

IP(10ROER(IIP*00) .EQ.O.ANO.IOROERtllPtuJP) .EQ.O)  60  TO  337 
IF(IORUER(ItP*00).'6T.O.ANO*IOROER(IXP*OOP).GT.O)  60  TO  335 
IP(I0HUER(IIP*001 .EQ.O)  60  TO  33A 
IITEMPallP 
OOTEMPaJO 
ICQOEse 
60  TO  337 
33A  IITEMPallP 
OuTEMPcJOP 
lOCOEsO 
60  TO  337 

335  60  TO  (336v3A2i339)  tiCOOt 

336  lAalURUER(IIiJu) 
iHalOROERdlPtUO) 

ICalORUERdIPvOOP) 

IOalOROER(II«OOP> 

II«IIP-1 

60  TO  336 

337  IlPallP*! 

IFdIP.OT.NSTOP)  CALL  ERRUR(3*I1P«NST0P.0J) 

60  TO  333 
336  NEL«N£L*I 

IF(NEL.6T. 10001  CALL  ERR0R(2»NEL*1000»0) 

IX(NEL»X>»1A 

IX(N£L«2)aIQ 

IX(NEL«3i«IC 

IX(NEL«4>«10 

1X(N£L»5)>MTL 

ALPHA  (NED  aBLKANO  (VTL) 
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GO  TO  ( 345«-343«340t34l,344  )•  ICOOE 

339  CONTINUE 

C TRIANGULAR  ELEMENTS 
lAalOROERdlf  JJ) 
lB«IOROER{IlTEMP,JgTEMP) 

IC>IORUER(IIiJuP) 

lOoIC 

60  TO  338 

340  lAsIOROCHdItUJ} 

IBalOROERdlPiUO) 

IC»IORUEW(  ITEMP.JoTEMP) 

1D*IC 

I0CDES4 
60  TO  336 

341  lAalUHUERdlPtUO) 

IBaIORU£R(IIP*JJP) 

IC»lORUEKdITENP<UgTEMP) 

IDalC 

IIaIIP-1 
lOCOEvl 
60  TO  336 

342  lAalORUERdltJU) 

IBalOROERdiTEMPf  JvTEMP) 

IC>I0RUERdl»3UP) 

IO«IC 
GO  TO  238 

343  lAtlORUERdlijOP) 
iBalOROERdlTEMP.JvTEMP) 
lC»lURUERdIPfUUP) 

ID*1C 

lOCOEat) 

GO  TO  338 

344  1A«I0RU£R(IIP*00) 

IB«lORUEHdlP»OOP) 

IC»I0R0ER(IITEMP,J^TEMP> 

IO«IC 

Il»llP-l 
I0C0E«1 
GO  TO  338 

345  DO  360  1HTL»1*NMTL 
IFdl.uT.MATRILdHTLfSnQO  TO  360 
IFdI.GE.MATPlL(IMTLt3n60  TO  360 
lF«OU.LT.MATRlL(lNtL*AnQO  TO  360 
IFJJO.GE.MATPILdMTLiSnGO  TO  360 
MTL  • MATRILdMTLtl) 

GO  TO  380 
360  CONTINUE 

PRINT  2100«NEL«II»wJ 
HTL»1 

380  IX(N£L«5)«HTL 

ALPHA  (NED  mBLKANG  CPTD 
GO  TO  313 
400  CONTINUE 
NUMELmNEL 

480  IF(NUMNP.GT.IOOO)  6R1TE (6*2002) 


i 

J 


I 

$ 


■'{ 
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T. 

f 


C SET  NOUAU  POINT  TEI^PERATUHE  TO  REFERENCE  TEMPERATURE 


C**««*«*««**«**«*»***« «««*•••»**•****«•*•#•*«*«*«*«««*«*«•««**«**«««**** 


IF(NUMTC.NE.O)60  TC  550 
DO  500  Nsl«NLMNP 
500  T0N>»TKEF 
550  WRITE(biAOOO) 

4000  FOHMATUH  * lOHENO  POINTS  ) 
RETURN 


C 

1000  FORMAT  (5I5«F10.01 
1002  format  (4IS*3F10.0) 

1300  FORMAT  (IHOt  AOHNOOAL  POINT  DATA— NO.  OF  NODAL  POINTS*  .15) 

2000  FORMAT ( 104HI  I J NP  TYPE  H-ORDINATE  2-OROINA 

ITE  H LOAD  OR  DISPLACEMENT  2 LOAD  OR  DISPLACEMENT) 

2001  F0HMAT(2lS.ie.F12.1,F12.3.Fl4.3iE26.7«E24.7.IlS) 

2002  FORMAT  (35H  BAD  INPUT  • TOO  MANY  NODAL  POINTS) 

2100  FORMATdHO.  EHELEMENT  . I4.3X.23HMITH  (I*D)  COORDINATES (* I2» IM. « 13* 
1 31H)  HAS  BEEN  ASSIGNED  MATERIAL  1 ) 

END 
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C data—  tins  ROUTINE  READS  CAROS  AND  PUTS  THEM  ON  TAPE  NO, 9 

SUBROUTINE  DATA 
C0RM0N/0H65/  JSTOP 
OIPENSION  CARO(20) 

REMIND  9 
Nao 

IF  ( JSTOP  .EQ.  I ) 60  TO  400 
READ(5«1000)  CARD 
10  »IRlTE(«jil002) 

60  TO  iiOO 

100  READ(5«lOOOiENOa300)CARO 
200  «RITE(6tl003)CARD 
MHITE(9»1000)  CARD 
NaN^l 
60  TO  100 
300  ENUFILE  9 
KEalNO  9 
HRITE{b«l00A) 

JSTOP=l 

RETURN 

400  kRlTE (btlOOS) 

STOP 

1000  FORMAT(20A4I 

1002  FORMATUni  /l3X,2M10»aX,2H20.8X,2H30»8X.2H40t8Xt2H50* 

18X«2H60«8X«aF70»eXi2H80/SX«BOHl234S67890 12345678901234567890123436 
278901234567890123456789012345678901234567890) 

1003  FOHMAT(5X,20A4) 

1004  FORMATUH  , IIHENO  OF  DATA) 

1005  format UH  , lOhENO  OF  jOB) 

ENC 


MNINK—  THIS  IS  A UTILITY  ROUTINE 
SUBROUTINE  MMMX(I«Jl 

CQr'HON/TO/IMIN(100)*INAX(100)«UMIN(30) iJNAX (30) *RAXI •MAXJ.NMTLtNBC 
I *HlNItMINJ 

IF(U.LT,UMIN(H > UHIN(I)«J 
IF(0.6T.JMAX(I) ) UFAX(I)>J 
IF(I,LT.IMIN(J) ) IFIN(U)aI 
1F(I.6T.IMAX(J)|  IFAX(U)>I 

RETURN 

ENC 
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C NODt—  A SMALL  LTILITY  ROUTINE  WHICH  MAY  BE  NEEDED  BY  MESH  OR  POINTS 

FUNCTION  node (It J) 

C 

COPMON/T0/IM1N(100)  tIMAXUOO}  tOMINOO)  tUMAXOO)  tHAXItMAXOtNMTLtNBC 
I tMINI.MINJ 
C 

N0CE>0 

DO  lOi)  UJaltiw 
NSTAKTalMINPwU) 

NSTUPalMAX (Jw) 

DO  100  lUNSTAHTtNSTOP 
N00E«NU0E*1 

IF(JJ.tQ.sl.ANO.II«EQ.I)  RETURN 
100  CONTINUE 
C 

HETUHN 

ENC 
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1 

10 

2 

11 

3 

12 

50 

13 


5UB80UTINE  ERROR (NtMtLtK) 

GO  TO  (1*2*3) *^ 

' •*I5*'II» 

GO  TO  50 

?S5IaTp“1*.JeL  » . *15. ‘EXCEEDS  MAXNEU  - ».15) 

GO  TO  50 

JSJJaTM^^I'Up’-  ‘*15. ‘EXCEEDS  NSTOP  - •*I5*‘0d 

FORMAT (‘  ‘.‘EXECUTION  T£HMlNATED‘) 

STOP 

ENC 


‘*15) 


> ‘*15) 
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